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Abstract 


This  report  describes  the  outcome  of  an  ionospheric  D-region 
study  made  with  the  wave  interaction  experimental  technique  mainly  using 
the  all-digitalized  facilities  at  The  Pennsylvania  State  University. 
Theoretical  basis  and  the  experimental  evidence  of  the  presence  of 
atmospheric  gravity  waves  in  the  D-region  ionosphere  are  presented. 
Spectral  analysis  of  the  wave  interaction  data  indicated  that  the 
presence  of  waves  with  periods  ranging  from  tens  of  minutes  to  few 
hours  are  possible.  Critical  intercomparison  of  two  independently- 
developed  electron-density  profile  inversion  techniques  are  presented. 

The  results  showed  that  the  two  methods  led  to  similar  profile  indicating 
their  reliability  in  profile  inversion.  Finally,  temporal  variations  in 
D-region  electron  densities  are  presented  showing  diurnal  variations  as 
well  as  the  variability  between  normal  and  anomalous  conditions. 
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1 . Introduction 

This  report  describes  the  outcome  of  an  ionospheric  D-region  study  made 
with  the  wave  interaction  experimental  technique  for  the  Office  of  Naval 
Research  under  Contract  Number  N0Q014-75-C-0178,  Modification  No.  P00001. 

The  research  program  used  in  part  the  wave  interaction  facilities  of  the 
Ionosphere  Research  Laboratory  at  The  Pennsylvania  State  University  which  have 
reached  a sophisticated  level  of  data  procurement  and  processing.  Data  were 
digitally  recorded  on  magnetic  tape  covering  28  D-region  heights  simultaneously 
every  300  milliseconds;  continuous  measurements  were  frequently  made  over  a 
time  span  of  over  10  hours,  with  additional  daily  sampling  of  D-region  electron 
densities  at  shorter  time  intervals.  Matching  to  this  data  procurement  capability 
is  an  elaborate  and  diversified  data  processing  high-speed  computer  program 
tailored  to  individual  research  objectives  to  be  described  in  detail  in  sub- 
sequent sections. 

The  main  objectives  of  this  research  program,  which  have  been  carried  out 
relying  on  the  sophisticated  capability  of  the  wave  interaction  experiment,  can 
be  summarized  as  follows: 

(a)  Effects  of  Dynamic  Processes  of  Atmospheric  Wave  Motions  in  the  D-Region 

Ionosphere  - 

Various  perturbations  in  the  ionosphere  due  to  inf rasonic  waves , atmospheric 
gravity  waves,  tidal  oscillations,  etc.  can  account  for  the  variability  in  the 
phase  and  amplitude  of  radio  signals  in  communication  systems.  In  spite  of  the 
fact  that  these  phenomena  are  present  in  the  ionosphere  most  experimental 
measurements  neglected  these  effects  being  unable  to  detect  such  effects.  With 
the  high  resolution  in  time  and  height  the  wave  interaction  experiment  is  capable 
of  identifying  the  atmospheric  waves  in  the  D-region.  Such  knowledge  not  only  will 
add  to  the  understanding  of  the  dynamics  of  the  D-region,  but  also  to  better 
model  the  chemistry  since  it  is  needed  in  properly  removing  the  influence  of 
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atmospheric  waves  from  data  in  order  to  correctly  account  for  the  solar  control 
upon  the  D-region.  The  effects  of  gravity  waves  were  observed  using  Penn  State 
facilities  at  State  College  as  well  as  at  Arecibo  NAIC  Observatory,  Puerto 
Rico.  Theoretical  background  and  the  experimental  results  on  the  subject  are 
described  in  detail  in  Section  2 of  this  report. 

(b)  Intercomparison  of  Electron  Density  Synthesis  Methods  - 

A reliable  data  reduction  technique  is  an  essential  part  of  this  research 
program  since  it  is  desirable  to  produce  accurate  electron  density  profiles 
from  wave  interaction  data  profiles.  The  profile  inversion  (or  synthesis) 


method  which  IRL  has  been  using  is  based  on  an  overdeterministic  scheme  developed 
by  Dale  in  1973  and  has  been  proven  to  be  reliable.  However,  in  view  of  the 
fact  that  another  profile  inversion  method  developed  by  C.  H.  Shellman  at  NELC 
for  processing  VLF  data  happened  to  have  considerable  similarities  compared  to 
that  of  Dale,  it  was  considered  important  to  intercompare  the  two  methods 
jointly  with  NELC  to  determine  whether  the  methods  could  complement  each 


other  to  improve  the  reliability  of  the  synthesized  electron  densities.  This 
cooperative  effort  started  with  initial  exchange  of  technical  informations, 
studying  and  testing  of  each  others  method  independently,  and  wound  up  with  two 
parties  closely  working  together  with  Mr.  Shellman  visiting  IRL  for  two  weeks 
in  June  1976.  In  spite  of  the  fact  that  there  are  some  fundamental  differences 
in  their  approach , the  results  produced  by  the  two  methods  were  found  to  be 
strikingly  similar.  It  was  satisfying  to  find  the  two  methods  confirm  each 
other  in  their  reliability.  Details  of  these  methods  and  the  test  results 
are  described  in  Section  3 of  this  report. 

(c)  Continued  Routine  Operation  of  Wave  Interaction  Facilities  - 


n 


its  ultimate  capability  in  studying  the  D-region  ionosphere.  For  this  purpose 
continued  routine  operation  of  the  facility  was  conducted  to  determine  and 
to  show  its  capability  in  observing  and  resolving  temporal  variations  of  the 
D-region.  The  results  of  this  effort  are  described  in  Section  A of  this 
report. 
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2 . Effects  of  Dynamic  Processes  of  Atmospheric  Wave  Motions 

For  some  years  the  Wave  Interaction  Facility  of  the  Ionosphere 
Research  Laboratory  at  The  Pennsylvania  State  University  has  been  used  to 
study  the  D-region  of  the  ionosphere  through  the  measurement  of  electron 
densities.  During  this  period,  improvements  have  been  made  in  the  experi- 
mental facility,  data  collection  techniques,  and  data  processing  methods. 

As  the  resolving  power  of  the  whole  system  has  increased,  it  has  become 
apparent  that  certain  fluctuations  in  the  data,  which  might  once  have  been 
considered  undesirable  noise,  are  caused  by  dynamic  processes  in  the 
D-region.  Those  processes,  and  how  they  may  be  studied  using  wave  inter- 
action are  the  subject  of  this  section. 

Specifically,  this  section  shows  how  the  wave  interaction  experiment 
may  be  used  in  the  study  of  gravity,  or  buoyancy,  waves  in  the  D-region  . 
These  waves,  believed  to  be  generated  near  the  surface  of  the  Earth, 
carry  energy  into  the  upper  atmosphere,  and,  surprisingly,  what  is  considered 
a small  perturbation  in  the  troposphere,  will  grow  and  have  a significant 
effect  on  the  dynamics  of  the  upper  atmosphere.  This  characteristic  and 
others  are  discussed  in  subsection  2.1,  which  gives  a brief  discussion  of 


the  theory  of  atmospheric  waves.  This  subsection  also  considers  the  means 
by  which  atmospheric  fluctuations  result  in  detectable  electron  density 
variations  as  measured  by  the  wave  interaction  experiment. 

The  second  subsection  2.2  describes  the  methods  by  which  the  data  is 
collected  and  processed.  It  is  divided  into  three  parts:  2.2.1  discusses 

appropriateness  of  methods,  2.2.2  discusses  Fourier  analysis  as  a tool 
for  data  analysis,  and  2.2.3  describes  the  interactive  program  used  for 
data  analysis  and  display. 


5 


The  third  subsection  2.3  presents  preliminary  results  found 
by  the  means  described  in  the  second  subsection  and  relates  them  to  the 
theory  discussed  in  the  first.  The  need  for  further  additional  techniques 
to  reveal  more  of  the  information  held  by  the  data  is  demonstrated  and 
some  on-going  work  in  this  area  is  described.  Finally,  the  knowledge 
gained  in  this  research  is  used  in  a discussion  of  improvements  which 
could  be  made  in  the  facility  to  obtain  a more  complete  knowledge  of  the 
dynamics  of  the  D-region. 

2.1.1  Theory  of  Low  Frequency  Atmospheric  Propagation 

The  theory  of  atmospheric  waves  is  an  application  of  fluid 
mechanics,  which  has  its  basis  in  the  derivation  of  microscopic  properties 
of  a medium  which  should  be  treated  statistically.  The  atmosphere  is 
composed  of  a very  large  number  of  molecules  so  that  for  all  practical 
purposes  the  physics  of  individual  particle  motion  leads  to  a statistical 
descriptions  of  overall  behavior  which  is  so  nearly  accurate  that  the 
statistical  base  may  be  dropped  leaving  a set  of  governing  laws  which  are 
deterministic  in  nature. 

Molecular  motion  may  be  classed  according  to  the  degree  of 
organization  or  coherence.  Wave  motion  in  the  atmosphere  is  a very 
organized  kind  of  motion  in  that  certain  vector  quantities  describing 
wave  motion  may  be  predictable  over  hundreds,  or  thousands  of  kilometers. 
Heat,  on  the  other  hand,  is  a very  disorganized  form  of  atmospheric 
motion  since  it  is  random  motion  of  the  individual  molecules  and  may  be 
described  on  a macroscopic  level  by  a single  scalar  variable,  temperature. 
Then  turbulence  is  a form  of  motion  which  is  organized  to  a lesser  degree 
than  wave  motion  but  more  so  than  heat,  since  it  may  be  described  as  the 
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more  or  less  random  motion  of  groups  of  molecules. 

One  reason  why  atmospheric  waves  are  of  interest  is  that  they 
carry  energy  from  one  place  to  another,  and  they  may  do  it  much  more 
rapidly  than  random  thermal  motion  (heat).  Energy  may  be  also  trans- 
ferred from  an  organized  form  to  a less  organized  one,  that  is,  waves  may 
generate  turbulence,  or  heat,  perhaps  in  a region  far  from  where  the  wave 
was  generated.  Since  the  waves  are  represented  by  variables  of  average 
motions  rather  than  individual  motion,  this  transfer  of  energy  is 
represented  by  a dissipation  term  in  the  equation  of  motion  much  as  it 
would  be  in  the  equations  for  a sliding  mass.  In  general,  dissipation 
is  very  difficult  to  model  mathematically  and  it  is  usually  ignored  if 
at  all  possible. 

Mathematical  modeling  of  the  macroscopic  properties  of  the 
atmosphere  leads  to  a set  of  differential  equations  which  in  general 
have  wave-like  solutions.  The  variables  and  parameters  in  the  equations 
may  be  divided  into  three  classes:  first,  those  which  are  directly 

derived  from  the  microscopic  properties  of  the  gas:  pressure,  density, 

heat  capacity  and  thermal  conductivity.  Second,  those  which  describe 
external  forces  acting  upon  the  gas  as  a whole,  such  as  gravity  and  the 
Coriolis  force.  Third,  those  variables  which  directly  describe  the 
response  of  the  medium.  This  is  usually  taken  to  be  a velocity  vector 
which  gives  the  macroscopic  or  average  velocity  of  the  medium  as  a function 
of  three  position  coordinates  and  time.  It  is  also  possible  to  use  a 
position  vector  which  describes  the  location  of  a particular  bit  of  air  as  a 
function  of  time  referred  to  an  initial  location. 
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The  mathematical  description  of  atmospheric  waves  involves  the 
use  of  three  equations:  motion,  continuity  and  state.  Certain  intuitive 

arguments  leading  to  these  equations  will  be  discussed  here;  for  more 
rigorous  derivations  refer  to  Gossard  and  Hooke  (1975).  The  approach 
adopted  here  is  known  as  the  Eulerian  in  which  the  motion  of  the  air 
is  described  by  its  velocity  as  indicated  above.  An  equivalent  description 
can  be  formulated  in  terms  of  the  position  vector;  this  is  known  as  the 
Lagrangian  system. 

The  equation  of  motion  is  a variation  on  the  usual  force  equals 
mass  times  acceleration  equation.  The  mass  is  replaced  by  the  density, 

P,  and  obviously  the  force  must  be  replaced  by  the  appropriate  quantity. 

The  force  per  unit  area  on  a wall  in  the  atmosphere  is  given  by  the 
difference  between  the  pressures  on  the  two  sides  or  by  the  gradient  of 
the  pressure,  p. 

The  two  external  forces  on  the  atmosphere  are  the  vertically  directed 
force  of  gravity  and  the  Coriolis  acceleration.  Actually  the  rotation 
of  the  earth  introduces  two  effects;  the  first  is  a latitude  dependent 
upward  acceleration  which  may  be  allowed  for  by  adopting  a modified  value 
for  g,  the  acceleration  due  to  gravity.  The  second  effect  is  a force 
directed  perpendicular  to  both  the  direction  of  motion  and  the  angular 
momentum  vector  of  the  earth.  The  latter  effect  may  be  ignored  for  the 
waves  with  periods  much  shorter  than  a day  to  be  considered  here. 

The  equation  of  motion  then  becomes: 

pg  - grad  p 2.1 
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For  static  conditions,  U = 0,  the  equation  becomes 


3p 

P080  " 3z 


where  the  subscript  zero  indicates  static  unperturbed  conditions,  and  z 

is  the  vertical  coordinate.  Elementary  gas  lows  show  that  the  density  and 

pressure  are  proportional  for  steady  conditions.  The  relationship  may  be 
2 

stated  as  Pq  = C /yPq>  where  C is  the  speed  of  sound,  about  300  m/sec,  and 

y is  the  ratio  of  constant  pressure  specific  heat  to  constant  volume 

specific  heat.  Substitutions  and  solution  of  the  first  order  differential 

2 

equation  shows  that  both  Pq  and  Pq  are  proportional  to  exp  (-ygz/C  ). 

2 

C /yg  is  defined  as  H,  known  as  the  'scale  height'. 

The  exact  differential  equations  which  results  from  the  equations 
being  considered  is  difficult  to  solve.  Thus  perturbations  causing 
wave  motions  are  assumed  small,  and  the  perturbed  value  is  set  equal  to 
the  static  value  plus  a perturbation.  This  is  substituted  into  the 
appropriate  equation  and  only  terms  up  to  the  first  order  are  kept.  Thus 
Equation  2.1  becomes 


pg  - grad  p 


where  p and  p are  the  perturbed  values. 

The  continuity  equation  is  an  expression  of  the  conservation  of  matter 
and  thus  relates  the  partial  derivitive  with  respect  to  time  of  the  density 


to  the  causes  of  density  variation.  First,  if  the  density  function  has 
spatial  changes,  that  is,  a non-zero  gradient,  then  velocity  of  the 


medium  will  cause  a net  transport  of  matter  to  or  from  any  given  location 
in  space.  Also  if  the  velocity  at  a particular  point  has  divergence  then 
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there  is  a net  accumulation  or  depletion  of  matter  at  this  point.  The 
equation,  including  only  first  order  perturbation  as  described  above  is: 


— + U • grad  pQ  + pQ  div  U = 0 


In  the  derivation  above  of  the  relationship  between  unperturbed 
density  and  pressure,  thermal  equilibrium  was  assumed.  For  wave-like 
perturbations,  the  pressure  and  density  changes  generally  flow  faster 
than  heat,  requiring  a relationship  which  is  obtained  when  there  is  no 
heat  flow  (adiabatic).  Some  simple  thermodynamics  (Gossard  and  Hooke, 
1975)  yields  the  relationship  . For  the  total  derivative, 

we  may  substitute  the  partial  derivative  plus  the  term  resulting  from  a 
non-zero . gradient . Dropping  the  appropriate  terms,  to  give  first  order 
perturbation  only,  gives 


grad  p 


? 

o-c  L 


3p/3t  + U • grad  p 


By  algebraic  manipulation  of  the  three  basic  Equations  2.3  through 
2.5,  they  may  be  modified  to  give  a form  for  which  a wave  solution  may 
be  found.  It  is  helpful,  however,  to  first  look  at  some  of  the  simple 
motions  possible  in  such  a medium,  and  it  is  especially  useful  to  look  at 
what  happens  if  some  of  the  terms  are  ignored. 

Consider  a medium  which  has  the  unperturbed  density  and  pressure  varia- 
tions of  the  atmosphere,  but  is  composed  of  an  incompressible  fluid.  Suppose 
a volume  element  of  this  fluid  is  given  a downward  initial  velocity.  This 
element,  with  density  , encounters  an  upward  restoring  force 
approximately  proportional  to  its  displacement  since  it  moves  into  a region 
of  increasing  density.  This  causes  the  element  to  move  back  up  and 
experiences  oscillation  about  its  initial  position. 
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To  find  the  frequency  of  this  oscillation,  the  appropriate 
relationship  must  be  substituted  into  the  equation  of  motion.  In  an 
incompressible  fluid,  the  divergence  of  the  velocity  is  zero,  and  the 
gradient  of  the  pressure  is  due  only  to  two  causes,  the  gravitational 
force  and  the  movement  of  the  volume  element  into  regions  of  different 
density.  The  first  of  the  two  must  cancel  the  pg”  term  in  Equation  2.3; 
this  is  how  the  exponential  height  variation  was  derived.  Thus  we  are 
left  with 

Po  It  = 'grad  P 2-6 

» 

where  p is  the  pressure  perturbation  due  to  the  movement  of  the  volume 
element.  If  the  element  always  has  the  same  density  as  the  medium  around 
it,  equilibrium  would  never  be  disturbed,  and  there  would  be  no  pressure 
perturbation.  If  the  density  were  zero,  the  pressure  on  the  element 
would  be  due  to  the  full  unperturbed  pressure  gradient.  In  general,  it 
will  be  given  by  the  fractional  density  change  times  the  pressure  gradient, 
that  is,  the  restoring  pressure  will  increase  as  the  element  moves  into 
regions  of  increasingly  different  density.  Thus  for  small  displacements. 
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grad  p'  - - - g A exp  ("Z/H)  z = - A exp  (~z/H)  . z 


when  z is  the  coordinate  of  vertical  displacement  such  that 


3z 

3t 


Then,  substituting  into  2.6  gives 

aZT* 

A exp(-z/H)  — y-  = A exp  (-Z/H)  z' 

3t  C 


2.7 


2 2 

3 z 1 _ g Yz1 

2 “2 
3t  C 


2.8 


A sinusoidal  solution  z = B sin  wt  gives 


w = g y 2.9 

C 

Allowing  compressibility  in  the  medium  modifies  the  frequency 
of  oscillation.  We  assume  that  the  process  is  adiabatic,  that  is,  no 
heat  is  lost  from  the  element  of  air  as  it  is  compressed.  This  simplifies 
the  analysis  and  is  a good  approximation.  If  the  process  were  isothermal, 
that  is,  heat  flowed  until  thermal  equilibrium  was  reached,  then  there 
would  be  no  oscillations.  We  can  see  this  by  imagining  an  element  of 
air  moved  downward  so  slowly  that  the  rise  in  temperature  from  compression 
is  negligible  as  heat  is  conducted  away.  Under  these  conditions  the 
product  of  pressure  and  volume  is  constant  and  thus  so  is  the  ratio  of 
pressure  and  density.  Since  the  energy  gained  from  compression  is 
conducted  away  as  heat  there  is  none  available  to  develop  a restoring 
force;  the  volume  element  of  air  is  always  identical  to  and  thus  inter- 
changeable with  its  ambient  elements. 

For  the  adiabatic  compression  case,  the  perturbed  pressure 
gradient  per  unit  displacement  will  be  given  by  the  ratio  of  density 


, •**“■’*  .VS*»  <> 


T*  rvyBK.-  ■ 


12 


gradient  to  density  times  the  unperturbed  pressure  gradient  as  in  the 
incompressible  case.  However,  the  density  gradient  must  now  be  modified 
by  the  adiabatic  compression  received  as  the  volume  element  moves  into 
regions  of  higher  pressure.  When  considering  adiabatic  changes  from  initial 
conditions  p,  p,  and  T the  relationship  which  holds  is 
_2  J 

dp  = C dp 

The  application  of  this  is  very  simple  since  it  is  only  necessary  to 
decrease  the  pressure  gradient  used  in  the  incompressible  case  by  an  amount 

= (1/C2)  ^ = (1/C2)  g A exp  (~Z/H)  = -^  A exp  (_Z/H)  2.10 

The  fractional  density  change  per  unit  displacement  is  then 


» 

i 


This  is  known  as  the  Vaisala-Brunt  frequency  (for  temperature  independent 
of  height)  and  it  appears  as  a parameter  in  most  equations  describing 


atmospheric  motion  when  the  density  and  pressure  gradients  caused  by  the 
force  of  gravity  are  considered. 

That  this  frequency  is  independent  of  pressure  and  density  may  seem 
remarkable.  It  should  be  recognized  that  any  resonance  may  be 
thought  of  as  the  interaction  of  a spring-like  restoring  force  and  mass 
or  inertia.  In  this  case,  the  physical  forces  governing  the  variations 


3 


in  one,  also  control  the  other,  leading  to  resonance  which  is  independent 
of  height  to  the  extent  that  the  temperature  (hence  C)  and  g remain 
constant . 

The  oscillation  described  here  represents  the  free  response  of  a 
small  bit  atmosphere.  It  is  not  a form  of  wave  propagation,  but  might 
be  a part  of  a process  in  which  wave  propagation  occurs.  If  there 
exists  some  source  of  random  atmospheric  motion  in  the  vertical  direction, 
the  maximum  velocity  and  displacement  in  the  response  would  be  expected 
at  ui  if  the  energy  of  the  source  is  approximately  uniformly  distributed 

D 

with  frequency  in  that  range.  However,  whether  or  not  wave  propagation 
occurs  at  <d„,  or  at  other  frequencies  may  be  most  easily  determined  by 

D 

solving  the  three  equations  for  wave  solutions. 

The  equation  of  motion  is  a vector  equation  and  thus  represents 
three  scalar  equations.  This  may  be  simplified  to  two  equations  without 
loss  of  generality  by  assuming  the  wave  propagation  in  the  horizontal 
occurs  only  in  one  direction,  say  the  x.  The  other  two  equations  are 
scalar,  giving  a total  of  four.  Some  simplification  may  be  gained  by 
assuming  zero  divergence  (incompressibility)  in  the  equation  of  continuity. 
The  loss  of  the  spring-like  restoring  force  provided  by  compressibility 
eliminates  from  the  solution  all  acoustic  waves,  but  not  the  so-called 
gravity  waves.  This  would  be  expected  from  the  above  analysis  in  which 
the  buoyancy  force  was  shown  to  be  modified  by,  but  not  dependent  upon, 
compressibility.  A solution  for  gravity  waves  in  the  absence  of 
compressibility  is  given  in  Gossard  and  Hooke  (1975). 
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For  the  scalar  equations  derived  from  the  three  original  equations, 

Hines  (1960)  assumed  a wave-like  solution  for  the  pressure  p,  density  P,  hori- 
zontal velocity  U and  vertical  velocity  U . Note  that  p and  P are  perturbed 
X z 

quantities,  while  p^  and  Pq  are  the  corresponding  steady  state  values,  and  the 
quantities  (p  - p^)  and  (p  - p^)  are  perturbations.  The  assumed  solution  is 
P * (p  - Pq)/p0  = R ' (P  ~ Pg)/po  = Ux/X  = u2 /z  = A exP  i(<*>t-  k^x  - f^z) 

2.13 

P,  R,  X and  Z are  complex,  and  give  the  relationships  among  the  several 
variables,  while  A is  an  arbitrary  amplitude,  and  are  the  horizontal 
and  vertical  complex  wave  numbers,  and  m is  the  angular  frequency  of  the  wave. 

Substitution  of  the  assumed  solution  into  the  four  equations  yields  four 

simultaneous  linear  homogeneous  equations  of  four  variables.  It  may  then 

be  shown  that  the  solution  is  valid  only  when  the  quantities  w,  K^,  and 

are  related  in  a manner  which  may  be  described  mathematically  by  a quadratic 
2 

equation  of  w . This  equation  is  known  as  the  dispersion  relationship; 
properly  interpreted,  it  gives  information  about  the  possible  directions  of 
propagation,  and  the  phase  and  group  velocities.  Given  that  the  requirements 
of  the  dispersion  relationship  are  met,  the  solution  yields  P,  R,  X,  and  Z 
in  terms  of  to  y K^,  K^,  and  the  parameters  in  the  three  originial  equations. 

These  equations,  which  give  the  relative  magnitude  and  phase  of  the  variables 
are  called  polarization  relations. 

The  absence  of  dissipation  terms  in  the  basic  equations,  and  the  constancy 
of  the  unperturbed  variables  in  the  horizontal  direction  make  the  following 
assumption  plausible:  the  wave  amplitude  suffers  no  decay  in  the  horizontal 

direction.  Thus  Hines  (1960)  assumed  that  Kx  is  real  and  found  that  the  dispersion 
relation  gives  two  possible  classes  of  waves.  One  of  these  is  that  Kz  is  purely 
imaginary  and  thus  energy  propagates  only  in  the  horizontal  direction  and  there 
is  no  phase  variation  in  the  vertical  direction  (surface  waves) . This  class 
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if 

was  not  of  interest  to  Hines  at  the  time  although  it  does  exist  in  the  ionosphere 

(Valverde  , 1958,  Thome,  1968).  The  other  mode  describes  those  waves  called 

internal  gravity  waves;  for  this  class  K * k + i/2H  where  k is  real, and 

z z z 

vertical  phase  propagation  is  allowed. 

This  implies  that  the  vertical  amplitude  is  growing  at  an  exponential 
rate.  Thus  a small  perturbation  at  the  surface  of  the  earth  may  grow  until 
at  high  altitudes  it  represents  a large  percentage  of  the  unperturbed  pressure. 
The  equations  presented  here  will  no  longer  adequately  describe  the  propagation 

for  two  reasons:  first,  the  linearizing  assumptions  by  which  the  higher  order 

[ j 

terms  were  dropped  are  no  longer  valid,  and  second,  dissipation  becomes  very 
important  under  these  conditions. 

Since  all  waves  in  this  class  have  the  exponential  growth  factor  the 
assumed  solution  may  be  modified  to  read 

A exp  (z/2H)  • exp  i(wt  - k^x  - kzz)  2-14 

where  kx  and  k^  are  both  real  for  this  class  of  solution.  Thus  the  dispersion 
relation  becomes 


I 

i 

I 

i 

l 

t 


2 222  22  7 72 

(of  - a>  ) a /C  - <o  (k  + k ) + <4  k = 0 2-15 

a x z B x 

where  ujg  has  already  been  defined  and  wa  = yg/2C  is  called  the  acoustic 
cuttoff  frequency.  This  equation  may  be  rewritten  in  the  form 

2 

7 7 7 7 7 

kZ  = (uf  - (o  )/C  + (-|  - 1)  k 2-16 

z a z x 

0) 

which  clearly  shows  that  either  u>  < ui„,  or  w > w since  both  k and  k are 

B a x z 

real  for  internal  waves.  Further  rearrangement  yields: 


p>*r*  • 
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2 ? 9 

(oj  - ui  ) / C“ 
a 


(a)2  - UJ2)  /c2 


If  u > a)  this  equation  describes  an  ellipse  centered  at  the  origin 
with  major  axis  a = 2 . ^ . As  u becomes  very  large,  the  major 


\ C2 (l-u)2 /w2) 
d a 


and  minor  axis  becomes  equal,  and  the  relationship  between  k and  k is  a 

x z 

circle,  as  expected  for  sound  waves.  When  o>  < co  < m , the  equation  does 

d a 

not  describe  any  real  figure,  but  when  a>  < ou„,  the  result  is  a hyperbola 

b 

centered  at  the  origin  with  the  slope  of  the  asymtotes  given  by 

, n 2,  2. -1/2 

+ (1-Wg/u  ) 


The  k^  versus  k^  relationship,  with  T = — as  a parameter,  is 

shown  in  Figure  2.1;  parameters  from  the  middle  D-region  are  used  (to  = 

.016,  ioa  = .0195).  It  should  be  remembered  that  this  represents  the  case 

^ = k^  , Kz  = k + i/2H,  and  that  other  solutions  are  possible.  The  ellipses 

are  acoustic  waves;  the  hyperbolas  are  gravity,  or  buoyancy,  waves.  Note 

that  the  direction  of  phase  propagation  is  given  by  a line  drawn  from  the 

origin  to  a point  on  one  of  the  curves,  and  that  for  long  period  buoyancy 

waves  this  direction  is  almost  constant  for  anv  k . This  implies  that 

z 

the  direction  of  propagation  is  almost  entirely  determined  by  the  frequency, 
and  also  it  is  clear  from  the  figure  that  the  horizontal  wave  number  is 
almost  fixed,  while  the  vertical  may  vary. 

In  order  to  look  for  other  solutions,  it  is  convenient  to  begin 

with  Hines  (1969)  original  dispersion  equation. 

42229  2o  2 

u»  - u C (K~  - K2)  + (y-1)  g2  K2  + iy gw  K = 0.  2.18 

X Z X Z 
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1/2 

Substituting  in  the  relationships  01  = yg/2C,  and  o>  = (y-1)  g/C 


gives 


22  222  222  2 4 

(Cuv,-u)C)K“-cj  C K + i2C<o  K = -oj 
n x z a z 


2.19 


Consider  the  case  (real),  and  K ^ = im^;  substituting  this  into  the 

above  equation  gives 

,„2  2 22,  . 2 . 2_2  2 2 4 

(C  io_  - to  C ) k + co  C m - 2Cio  to  m = -to 
B x z a z 

Completing  the  square  and  dividing  through  puts  this  equation  in  conic 
form: 


(m  -to  /CY 
z a 


(to2-to2)/C2 

a 


2 2 2 
(« -(o  )/<r 

a 


(l-^/to/) 


Note  that  for  large  to  this  is  a hyperbola,  and  for  small  to  it  is  an  ellipse, 

just  the  opposite  as  in  the  last  case.  This  case  represents  no  vertical 

phase  propagation  and  the  waves  are  called  Lamp  waves,  and  the  k^  - k^ 

relationship  is  shown  in  Figure  2.2. 

Consider  next  the  possibility  of  being  imaginary,  say  = im^. 

Then,  if  K = k -i/2H  as  for  internal  waves  the  onlv  difference  between  this 
z z J 

2 

case  and  that  of  the  internal  waves  is  that  the  sign  in  front  of  the  m 

x 

term  will  be  changed.  There  will  be  no  solution  for  <o  < to  since  there  can 

2 2 B 

be  no  conic  of  the  form  - = 1.  On  the  other  hand  <o  > co  will  give 

a b 

2 

X V 

hyperbola  of  the  form  - — + -“j  = !•  There  will  be  vertical  hyperbola 

a b 

rather  than  horizontal  as  before,  and  the  slopes  of  the  asymototes  will  be 
2 

given  by  + (1  - ^B) . This  relationship  is  shown  in  Figure  2.3. 

(o2  2 2 

Finally,  consider  K complex,  say  K = L + i M . Since  K =L  +2iL  M 

X XXX  XX  XI 

2. 

- M no  relationship  between  K and  K with  onlv  real  coefficients  mav  be 
X X z 

established.  Thus  there  can  be  no  waves  with  Kx  complex. 
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It  is  convenient  to  think  of  the  medium  having  an  index  of 

refraction  of  one  if  ui/k  = C.  In  general,  to/k  = C n,  where  n is  the  index 

k C 

of  refraction  and  n = — 1 — , n = k C/ou.  Figure  2.4  shows  the  gravity 
x to  z z 

wave  and  acoustic  wave  relationship  for  n^-n^. 

For  the  internal  gravity  waves  and  acoustic  waves,  the  direction 
of  energy  flow  is  not  readily  apparent.  From  the  wave  propagation  theory  it 
is  known  that  it  must  be  in  the  direction  of  the  group  velocity,  which  has 

r “j  -g 

the  components  i 3u)/3k  , and  i . Thus  the  direction  may  be 

L XJ  x L z_  x 

roughly  determined  by  looking  at  the  k^-k^  curves  of  constant  T.  As  Hines 
(1960)  shows,  this  is  equivalent  to  the  direction  of  a normal  drawn  to  the 
curve  away  from  the  origin. 

The  effect  of  dissipation  will  not  be  discussed  here  except  to 
say  that  higher  frequency  waves  are  effected  the  most  by  viscosity.  Thus 
very  short  wavelengths  will  be  absorbed  before  reaching  the  D- region.  On 
the  other  hand,  waves  with  long  vertical  wavelength  (low  frequency)  will 


be  reflected  by  temperature  gradients  unless  the  horizontal  wavelengths 
are  also  quite  long.  As  will  be  shown  later,  the  wave  interaction  experi- 


ment is  not  sensitive  to  the  vertical  structure  of  a wave  unless  the 
horizontal  structure  is  large. 


2.1.2  Description  of  Experimental  and  Data-Collecting  Facility 

The  experimental  and  data-collecting  facility  used  for  the  wave  interaction 
experiment  has  been  described  in  some  detail  by  Swanson  (1976).  The  description 
here  will  be  brief,  covering  the  basic  method  by  means  of  which  simultaneous 
data  is  generated  over  a height  range  covering  the  entire  D-region.  Sampling 
rate  and  data  accuracy  are  also  discussed. 
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Figure  2.5  shows  the  time  versus  height  relationship  for  a 
conceptual  version  of  the  experiment  involving  two  very  narrow  transmitted 
pulses,  the  W(wanted)  and  the  D (disturbing).  The  W is  relatively  low 
power  (5  kw  peak,  2.4  MHz)  and  is  transmitted  at  a rate  of  about  300/sec, 
while  the  D is  about  11  MW  ERP  (4.55  MHz)  and  is  transmitted  at  one-half 
the  rate  of  W pulse.  The  W pulse  reflects  from  the  E-layer  and  meets  the 
up-going  D pulse  at  a height  tr  which  depends  upon  the  relative  time  delay 
between  the  two  pulses.  The  set  of  pulses  (W^)  is  a set  of  successively 
delayed  W pulses,  which  meet  the  D pulse  at  successively  higher  tr  . 

The  experiment  may  be  thought  of  as  transmitting  the  set  (W^) , 
followed  by  the  D pulse,  and  then  receiving  the  series  of  echos  from  the 
set  (W^) . The  amplitude  of  each  of  the  W pulses  is  measured,  and  then  the 
sequence  is  begun  again,  only  this  time  after  transmitting  the  set  (W^), 
no  D pulse  is  transmitted.  The  amplitudes  are  measured  for  this  undisturbed 
cycle  and  compared  with  those  measured  during  the  last  cycle.  Cycles 
continue  to  alternate  in  this  manner,  and  the  comparison  between  the  two 
cycles  is  made  in  a device  known  as  the  wave  interaction  detector,  which 
gives  an  average  value  of  amplitude  differences  between  cycles.  Data  is 
taken  at  28  heights,  requiring  28  detectors,  followed  by  some  method  of 
storing  the  data. 

The  actual  experiment  differs  somewhat  from  the  above  simplified 
description.  Let  each  member  of  (W^)  be  one  kilometer  wide  with  one  kilo- 
meter spacing,  thus  creating  a single  wide  rf  pulse.  This  is  what  is 
actually  transmitted  and  known  as  the  W-pulse.  The  received  echo  is  divided 
into  the  proper  sections  in  the  receiving  equipment.  In  practice,  h^  is 
determined  from  the  delay  between  the  D pulse  and  the  leading  edge  of  the 
R pulse,  which  is  used  to  gate  the  returning  echo.  The  W pulse  is  made 
somewhat  wider  than  the  absolute  minimum  (which  is  the  width  of  the  R pulse) , 
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thus  allowing  some  leeway  for  variations  in  reflection  height. 

There  are  two  more  ways  in  which  the  practical  experiment  varies  from  the 
ideal.  First,  the  D pulse  is  finite  in  width  and  is  not  rectangular  in  shape, 
but  cosine  squared.  This  complicates  the  mathematics  of  data  reduction. 

Second,  the  received  echo  is  divided  into  overlapping  sections  wider  than  1 km, 
rather  than  distinct  1 km  sections.  This  results  in  some  smoothing  and 
eliminates  saturation  in  the  wave  interaction  detectors,  but  also  complicates 
data  analysis.  However,  these  variations  are  all  accounted  for  in  the  wave 
interaction  theory. 

The  data  to  be  recorded  consists  of  28  wave  interaction  detector  outputs, 
an  AGC  output  (Automatic  Gain  Control;  gives  the  amplitude  of  the  received  echo), 
and  a signalderived from  the  D-transmitter  monitoring  circuitry  which  gives  the 
transmitted  power.  The  30  signals  are  sequentially  sampled,  digitilized, 
serialized,  and  recorded  on  magnetic  tape  as  a first  step  in  the 
analysis  process.  The  sampling  rate  is  about  3 samples  per  channel  per  second, 
while  the  data  is  filtered  at  20  db /decade  with  a time  constant  of  about  10 
seconds.  Thus  the  sampling  rate  is  sufficient  to  obtain  highly  accurate  data  since 
aliasing  errors  will  be  very  small  indeed. 

The  data  is  digitized  to  a resolution  of  cne  part  per  thousand  full-scale. 

Since  data  values  generally  run  no  more  than  half  scale,  the  resolution  must  be 
decreased  by  this  factor.  On  the  other  hand,  means  obtained  by  summing  many 
points  will  have  resolution  increased  by  approximately  the  number  of  points 
summed.  This  will  be  true  as  long  as  the  sampled  signal  is  sufficiently  varying 

in  nature  so  that  the  least  significant  digit  changes  frequently  in  the  shortest 
observed  period.  This  condition  has  been  observed  to  hold. 
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2.1.3  The  Wave  Interaction  Experiment  for  Detecting  Atmospheric  Waves 

The  experimental  system  described  in  the  last  subsection  allows  continuous 
recording  of  data  which  is  related  to  electron  densities,  and  their  fluctuations. 
It  is  the  purpose  of  this  subsection  to  show  that  this  data  may  be  used  to 
determine  the  existence  of  atmospheric  waves  in  the  D-region.  There  are  two 
parts  to  this  problem:  first,  the  effect  of  the  passage  of  atmospheric  waves 

upon  the  electron  density,  and  second,  the  detection  of  this  effect  by  means  of 
the  wave  interaction  experiment. 

As  explained  in  Gossard  and  Hooke  (1975),  electrons  move  with  the  neutral 
particles  in  the  D-region  since  the  atmospheric  density  is  great  enough  so  that 
the  region  is  collision  dominated.  This  gives  a tremendous  simplicity  to  any 
analysis  as  compared  to  the  plasma  physics  required  in  regions  above,  but  also 
adds  one  significant  complicating  factor  for  longer  period  waves.  The  time 
constants  for  generation  and  decay  of  electrons  in  the  D-region  are  such  that 
longer  period  waves  may  modify  the  photo-chemical  reactions,  having  an  effect 
upon  electron  densities  which  is  unknown  at  present.  This  problem  is  quite 
complicated  in  theoretical  nature,  and  will  be  ignored  in  this  report. 

In  principle,  electron  density  fluctuation  could  be  detected  from  wave 
interaction  data  by  reducing  a sufficient  time  period  of  data  to  electron 
densities  and  then  searching  for  fluctuations.  However,  for  preliminary 
analysis,  the  data  values  themselves  serve  the  purpose  adequately  since  they  are 
approximately  proportional  to  electron  densities  with  the  restriction  that  the 
data  in  its  raw  form  does  somewhat  limit  height  resolution  because  of  the  finite- 
pulse  widths  and  nature  of  the  experiment.  However,  the  height  resolution  :s 
not  irretrievably  lost,  and  methods  which  allow  its  recoverv  without  cne  ;s * 
a complicated  synthesis  program  are  presently  under  development. 
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Although  the  wave  interaction  data  are  approximately  proportional 
to  electron  density,  it  is  also  dependent  upon  the  electron -neutral  collision 
frequency  which  depends  upon  the  neutral  pressure.  Atmospheric  waves  cause 
perturbations  in  the  pressure  and  thus  will  cause  variations  in  the  data 
values  through  modification  of  the  collision  frequency. 

The  approximate  significance  of  this  effect  can  more  readily  be 
evaluated  using  Appleton-Hartee  theory.  Equation  2-36a  of  Weisbrod  (1964) 
gives  a quantity  showing  the  effect  of  collision  frequency  variations  at  a 
particular  height 

fz  (1+Y)2-Z2] 

T = — ~ - — 2-20 

hA  h2  f(l  + Y )2  + Z2!  2 

u L _J 

where  Z = — — 

0) 

v = electron  collision  frequency 

to  = angular  wave  frequency 

A = constant  virtually  independent  of  height 

h = height 

Y = — cos  0 

L (o 

w = angular  gyrofrequency 
n 

0 = angle  between  the  direction  of  propagation  and  the  magnetic 
field  (0  **  19°  at  State  College). 

To  evaluate  this  effect,  and  compare  it  to  that  of  electron  density 
variations,  we  must  define  an  experimental  sensitivity.  For  example,  if  the 
wave  interaction  coefficient  T(h^)  is  assumed  to  be  linearly  proportional  to 
electron  densities  since  it  is  approximately  so  as  mentioned  earlier,  it 
can  be  given  by  T*  (h)  = K(h)  Ng  (h) . Then  an  acceptable  form  of  definition 
of  relative  sensitivity  to  electron  density  variations  is 
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giving  a direct  linear  dependency.  If  similar  form  of  definition  is  adopted 
for  the  sensitivity  due  to  collision  frequency 
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based  on  equation  2-20  where  D ® (1  ± 


This  turns  out  to  be  not  as  simple  as  Sjj  and  it  needs  to  be  calculated 


as  a function  of  height  for  typical  values  of  v.  The  results  are  shown  in 


Figure  2. 6 using  an  examplary  value  of  D = .61,  showing  that  S z is  reasonably 


constant  down  to  65  kilometers.  Defining  ST  * + SN  gives  the  total  relative 


sensitivity,  which  remains  reasonably  useful  to  somewhat  under  65  kilometers. 


At  lower  altitudes,  the  extremely  large  variation  of  S^,  with  h would  make 


the  results  very  difficult  to  interpret.  As  the  quasi-longitudinal  approxima- 
tion is  in  effect  here,  the  Sen-Wyller  equations  would  give  very  similar 
results . 

Another  factor  which  requires  careful  consideration  is  the  fact  that  the 
wave  interaction  experiment  is  affected  by  two  classes  of  noise:  first, 

interfering  r.f.  signals  on  the  wanted  frequency  channel,  and  second,  variations 
in  wanted  echo  strength  due  to  causes  other  than  wave  interaction.  Both 
of  these  sources  can  have  random  components  which  have  to  be  distinguished 
from  variations  in  alternate  pulse  area  due  to  wave  interaction. 
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The  first  class  consists  of  natural  atmospheric  noise  such  as 
lightening,  noise  spikes  generated  by  power  lines,  and  occassionally , man-made 
r.f. signals.  This  class  of  noise  is  significant  only  when  the  echo  is  weak, 
and  thus  tends  to  be  somewhat  more  important  at  noon  than  in  the  morning 
or  afternoon.  Around  sunset,  however,  noise  may  increase  tremendously  due 
to  changing  propagation  conditions  and  is  significant  even  when  the  echo  is 
very  strong. 

The  second  class  is  mainly  the  result  of  ionospheric  variations  near 
the  reflection  height.  Although  such  noise  is  present  all  the  time,  it  is 
most  important  during  the  times  the  E-layer  is  forming  or  decaying,  while  the 
echo  strength  tends  to  be  relatively  constant  around  noon.  Power  line  fluctua- 
tions and  power  supply  ripple  in  the  wanted  transmitter  contribute  no 
significant  noise  since  the  former  are  removed  by  the  AGC,  and  the  latter  are 
too  far  removed  from  the  frequency  of  coherent  detection  to  be  significant. 
Variations  in  D-transmitter  power  are  removed  from  the  data  by  recording  the 
actual  transmitted  power  and  then  normalizing  the  data  to  a standard  power. 

It  is  essential  to  know  if  variations  in  data  resulting  from  noise 
will  mask  variations  from  dynamic  behavior.  A very  simple  technique  has  been 
adopted  to  check  this,  which  involves  executing  the  experiment  with  the 
D-transmitter  turned  off.  The  W transmitter  is  operated  normally  and  both 
classes  of  noise  variations  will  be  present.  The  resulting  noise  data  is 
treated  just  as  if  it  were  real  data,  and  a power  spectrum  is  computed.  This 
is  compared  to  a power  spectrum  of  real  data,  giving  a signal  to  noise  ratio 
as  a function  of  wave  frequency.  In  this  manner,  the  useful  frequency  range 
of  the  experiment  may  be  determined.  This  data  will  be  presented  later  in  the 
report. 
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If  it  is  assumed  for  a moment  that  the  D-transmitter  uniformly  neats  a 
cone  of  the  ionosphere  with  half-angle  0,  as  shown  in  Figure  2.7(  then  a wave 
propagating  through  the  heated  region  with  wave  numbers  and  will  perturb 
the  electron  density  by  an  amount  proportional  to  the  amplitude  of  the  wave. 
Since  perturbations  of  data  value  at  a particular  height  represent  the  integra- 
ted effect  across  the  entire  arc,  the  relative  strength  of  the  perturbation 
divided  by  the  unperturbed  condition  is 


1 

P(z)  = / (N(h)  • A exp  j (K  x + K z - wt)  ) dx 

^ x z 

i 

/ N(h)  dx 

-1 

where  1 = Z tan  9,  and  the  arc  has  been  approximated  by  a line. 

This  reduces  to 

|P(z)  | = A • exp  j (K  z - wt)  (sinfk  z tan  3)/(K  z tan  9) 

which  shows  that  the  sensitivity  to  measuring  vertical  wave  structure  is 
affected  by  the  horizontal  structure,  the  altitude  and  the  size  of  the  heating 
cone.  Actually,  the  cone  is  not  abrupt,  and  the  mathematical  relationship 
governing  its  intensity  should  be  included  in  the  above  integral.  Despite 
this  approximation  and  the  approximation  made  in  the  integration,  it  is  clear 
that  small  horizontal  wavelength  (large  K^)  will  cause  the  loss  of  sensitivity 
regarding  the  vertical  structure.  For  very  small  K (large  X ), 

X X 

(sin|  z tan  9j)/(t<x  ztan  9)  approach  unity  and  there  is  no  loss  of  sensitivity. 

As  stated  earlier,  only  certain  waves  reach  the  D-region,  due  to  reflection 
or  dissipation.  Figure  6 in  Hines  (1960)  shows  that  waves  with  a vertical 
wavelength  of  three  or  more  kilometers  could  have  a wide  range  of  horizontal 
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wavelength  up  to  a few  hundred  kilometers  at  an  altitude  of  90  km.  On 
the  other  hand,  waves  of  vertical  wavelengths  of  approximately  up  to  70  km 
could  exist  at  90  km.  In  view  of  the  fact  that  the  experimental  capability 
in  detecting  wave  structures  is  ideally  around  vertical  wavelengths  of  a few 
kilometers  and  horizontal  wavelengths  over  few  tens  of  kilometers,  it  is 
considered  well  suited  for  studying  D-region  atmospheric  waves. 
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2.2.1  Methods  of  Data  Analysis 

The  first  step  in  data  analysis  is  to  determine  the  power  spectrum  of  data 
fluctuations  and  compare  it  to  the  spectrum  of  noise  fluctuations.  This  is 
relatively  simple  but  is  indispensable  in  evaluating  the  significance  of  the 
results.  The  second  and  third  parts  of  this  subsection  describe  the  analysis 
procedures  while  the  results  are  presented  in  2.3. 

The  previous  section  of  this  report  has  shown  how  the  wave  interaction 
experiment  has  the  potential  to  explore  the  vertical  dynamic  structure  of  the 
D-region,  with  some  inherent  limitations.  Demonstrating  the  existence  of  waves 
requires  considerable  analysis.  The  first  step  in  a thorough  analysis  is  to 
reduce  the  data  so  as  to  obtain  the  highest  possible  height  resolution.  The 
algorithm  for  this  process  is  nearly  complete,  but  will  not  be  described  here 
due  to  space  limitations.  The  major  limitation  in  the  process  is  that  a 
lower  boundary  condition  must  be  assumed  and  the  accuracy  of  the  results  does 
depend  upon  the  suitability  of  the  boundary  condition.  Although  data  is 
recorded  every  kilometer  over  much  of  the  height  range,  the  height  resolution 
is  estimated  as  two  to  three  kilometers  throughout  the  height  range. 

After  this  process,  the  data  may  be  thought  of  as  a two  dimensional  array, 
with  variables  of  height  and  time.  If  a few  strong  waves  present  at  discrete 
frequencies,  the  autocorrelation  function  with  height  as  a variable  would  show 
them  up.  However,  it  is  unlikely  that  such  a simple  situation  pertains,  since 
the  spectrum  is  likely  to  be  broad  and  continuous.  Cross  correlation  functions 
between  the  data  measured  at  different  heights  with  time  as  a variable  can 
enlighten  such  a spectrum.  However,  as  pointed  out  in  Gossard  and  Paulson  (1968) 
the  variations  in  phase  speed  will  tend  to  destroy  the  significance  of  the  cross- 
correlation in  the  three  receiver  experiment  for  detecting  horizontal  infrasound 
propagation.  Similarly,  the  variations  in  vertical  phase  speed  in  this  experiment 
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will  have  the  same  effect.  For  this  reason  Gossard  relied  on  the  so-called 
cospectrum,  which  may  be  loosely  interpreted  as  a cross-correlation  over  a 
limited  frequency  range.  Nevertheless,  data  analysis  must  proceed  with  extreme 
caution  in  a case  such  as  this. 


and  it  is  the  intent  here  to  show  that  the 


quantity  may  be  interpreted  in  such  a manner  as  to  yield  a quantitative  picture 


of  the  periodic  fluctuations  which  compose  the  signal  X(t).  Since  the  difficulties 


may  not  be  immediately  apparent,  a few  problems  will  be  described  here.  The 


approach  here  will  be  practical  in  nature;  the  theory  may  be  found  in  Blackman 


and  Tukey  (1959) 


The  most  easily  manageable  situation  would  seem  to  be  a signal  composed  of  a 


few  discrete  sinusoidal  waves 


This  function  is  shown  in  Figure  2.8,  and  its  power  spectrum  in  Figure  2.9 
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Although  three  peaks  are  clearly  distinguishable,  in  a real  situation  there 

would  be  some  question  as  to  whether  the  smaller  peaks  represent  real  data  or; 

as  in  this  example,  the  samller  lobes  of  three  overlapping  sine  functions 
s in  x 

( — - — ).  It  is  thus  desirable  to  modify  the  power  spectrum  so  that  a single 
period  present  in  the  data  will  yield  a single  peak,  even  if  some  resolution 
is  lost  in  the  process. 

Physical  data  does  not  usually  exist  without  a large  amplitude,  slow- 
varying  essentially  non-periodic  part,  which  we  shall  denote  as  a trend.  Such 
a trend  has  been  added  to  the  data  of  the  last  example,  and  the  result  is  shown 
in  Figure  2.10.  As  expected,  the  trend  has  a severe  effect  on  the  spectrum 
(Figure  2.11)  as  can  be  seen  in  the  figure. 

The  second  class  of  problems  has  to  do  with  the  integration.  First,  the 
function  is  not  known  continuously,  but  only  at  a set  of  discrete  points,  and 
second,  even  if  it  were,  the  integration  would  have  to  be  performed  numerically, 
and  thus  approximately,  since  no  analytic  representation  is  possible.  These 
two  limitations  complement  one  another,  and  an  error  bound  may  be  deduced 
readily. 

Numerical  integration  is  usually  carried  out  by  sampling  the  function  at 
evenly  spaced  points , finding  an  analytic  approximation  between  these  points  and 
adding  together  the  approximated  areas.  The  data  collection  process  yields 
evenly  spaced  samples  with  spacing  very  much  closer  than  the  periods  this  work 
is  concerned  with  here.  As  shown  earlier,  some  averaging  may  be  performed 
giving  a set  of  high  resolution,  equally  spaced  points,  with  spacing  much  wider 
than  the  original  data,  but  still  very  much  less  than  periods  of  interest. 

This  latter  set  of  points  is  the  one  used  for  Fourier  analysis. 

The  error  in  the  power  spectrum  may  now  be  related  to  the  error  inherent  in 
the  numerical  integration.  The  routine  used  here  employs  trapezoidal  integration, 
for  which  the  error  bound  is  given  by  *■ 
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e <_  (max  X (t)) 


where  T is  the  length  of  the  interval  of  integration,  m is  one  less  than  the 
number  of  points,  and  max  X(t)"is  the  maximum  value  of  the  second  derivative 
of  the  signal  in  the  interval  (Moursand  and  Davis,  1967). 

Assume  that  the  averaging  has  given  a signal  that  is  approximately 
bandlimited  to  gjo.  This  is  equivalent  to  saying  that  the  signal  changes  much 
more  slowly  than  a sine  wave  of  amplitude  approximately  equal  to  the  signal 
amplitude.  The  first  derivative  of  the  real  part  of  the  integrand  of  equation 
2-22  is  given  by 


X'  (t)  = - to  X(t)  sin  tot  + d cos  to(t) 

at 


In  looking  for  the  maximum,  the  second  term  will  be  negligible  and  may  be  dropped. 

2 

Using  the  same  approximation,  the  second  derivative  is  given  by  -to  X(t)cos  io(t). 


| - to^  X(t)  cos  to(t)  | | to^  max  X(t) 

It  is  convenient  to  set  T = 1,  and  so 


to  max  X (t) 


(f/m)  max  X(t) 


Thus  the  relative  error  to  the  maximum  of  X(t)  is  just  (f/m)  for  the  real  part 
of  the  integrand.  The  imaginary  part  will  have  the  same  error;  squaring  and 


adding  gives  an  error  in  the  power  spectrum  of 
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Note  that  a frequency  of  1 has  one  full  period  in  the  length  of  the  data  sample. 

As  long  as  we  have  at  least  four  points  in  the  period  of  the  highest  frequency 

4 

of  interest  the  spectrum  will  be  quite  accurate  since  2(f/m)  <_  .0078.  This  is 

just  the  error  inherent  in  the  mathematics;  the  errors  of  interpretation  due  to 
leakage  and  trends  described  earlier  must  still  be  dealt  with. 

Slowly  varying  trends  may  be  removed  by  fitting  a polynomial  to  the  data. 

For  the  data  in  this  report  a fourth  order  polynomial  fitted  by  the  least 
squared  error  method  was  chosen.  The  trend  varies  sufficiently  rapidly  that 
the  diurnal  variation  of  the  atmosphere  may  be  followed,  but  not  so  quickly 

that  useful  data  will  be  filtered  out.  The  spectrum  of  the  data  of  Figure  2.10 

is  shown  after  trend  removal  in  Figure  2.12.  Note  that  it  is  very  close  to  the 
spectrum  obtained  before  the  trend  was  inserted  (Figure  2.9). 

The  leakage  may  be  eliminated  by  a process  known  as  "Hanning"  which  consists 
of  replacing  each  coefficient  in  the  real  and  imaginary  part  of  the  transform 
by  a weighted  average  of  neighboring  coefficients.  The  theory  of  this 
operation  is  explained  in  Blackman  and  Tukey  (1959).  The  effect  is  shown 
in  Figure  2.13,  where  the  bar  indicates  the  range  of  significantly  weighted 
averaging.  Note  that  complete  elimination  of  the  minor  lobes  involves 
broadening  of  the  major  peak. 

2.2.3  The  Interactive  Program  for  Data  Analysis  and  Display 

The  Hybrid  Computer  Laboratory  of  the  Electrical  Engineering  Department 
of  The  Pennsylvania  State  University  operates  a PDP-10  computer  in  time-shared 
mode.  The  computer  has  a high  speed  CPU  and  a large  memory,  as  well  as  interfaces 
and  software  for  various  types  of  terminals,  display  screens  and  plotters. 

Used  in  conjunction  with  the  PSU  computation  center's  main  computer,  the 

I3M  370-168,  the  PDP-10  forms  a powerful  and  versatile  tool  for  data  analysis  and 

display. 
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A full  day's  operation  at  the  Scotia  experimental  facility  produces  about 
4 million  data  points  on  a seven  track  magnetic  tape.  The  370-168  is  used  for 
the  brute  force  averaging  into  sets  of  these  points  for  permanent  storage  in 
condensed  form  and  for  transfer  to  the  PDP-10  for  further  analysis.  This 
subsection  describes  the  interactive  program  used  in  this  analysis. 

There  are  several  basic  functions  which  the  program  must  perform.  These 
are  the  following:  (1)  the  reading  of  data  from  tape  or  disk,  (2)  the  processing 
of  the  data  as  described  in  the  last  subsection,  and  (3)  writing  out  or  displaying 
the  data,  both  the  final  result,  and  at  intermediate  points  to  aid  and  direct 
further  processing.  The  second  category  consists  of  various  types  of  averaging 
and  filtering  and  Fourier  analysis,  while  the  third  category  includes  numerical 
printouts,  display  screen  graphs,  and  plotter  graphs. 

The  strength  of  a time  shared  computer  system  lies  in  the  interaction  between 
the  operator  and  program  during  execution  time,  since  the  operator  is  able  to 
look  at  intermediate  results  and  make  decisions  which  are  beyond  the  capability 
of  the  machine.  However,  the  program  itself  must  be  designed  in  a manner 
which  increases  the  extent  of  this  interaction  by  making  the  most  information 
available  to  the  operator  in  the  most  convenient  manner. 

The  program  is  primarily  intended  to  be  operated  from  a Tektronix  4010 
graphics  display  terminal,  although  many  of  its  features  may  be  used  with  any 
terminal.  The  program  consists  of  a set  of  commands  which  may  be  executed  from 
the  terminal,  and  the  software  for  executing  these  commands.  These  commands 
each  perform  in  one  of  the  three  basic  categories  described  above.  When  execution 
of  a command  is  completed,  the  program  asks  for  another  command,  unless  instructed 
to  stop. 

A complete  description  of  the  program  would  be  too  long,  and  so  a few  examples 
will  suffice  here.  The  figures  in  the  last  subsection  were  produced  by  the 
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following  method:  first,  a command  denoted  DATA  was  executed,  computing  a 

test  function  and  loading  it  into  an  array  for  further  processing.  Using 
initial  information  regarding  the  frequency  range  of  the  function  certain  parameters 
are  set  using  the  command  SPAR.  If  at  any  time  during  further  processing  the 
initial  assumptions  change,  the  parameters  may  be  reset.  Then  the  command 
AVER  filter  out  unwanted  high  frequencies  and  reduces  the  number  of  points  in  a 
manner  set  by  SPAR,  to  enable  more  efficient  use  of  the  plotting  and  analysis 
routines.  The  function  is  then  plotted  on  the  graphics  terminal  to  check  for 
possible  errors,  and  finally  the  plot  is  made  on  the  HP-7202A.  In  a similar 
manner  commands  for  performing  trend  removal,  Fourier  analysis,  Hanning,  and 
so  on  are  performed  and  the  results  are  checked  in  the  graphics  terminal  and 
then  plotted. 

The  three-dimensional  plots  presented  in  the  next  subsection  were  produced 
by  a somewhat  more  elaborate  process.  First  the  data  is  read  and  trial  processing 
is  performed  and  checked  on  the  graphics  terminal.  When  the  proper  parameters 
and  procedures  have  been  decided  upon,  the  proper  command  sequence  is  stored 
and  executed  iteratively  on  successive  records.  The  plot  is  not  performed  on 
line,  but  is  stored  on  magnetic  tape  played  back  later  into  the  plotter  with 
the  assistance  of  IRL’s  NOVA  computer. 

2.3.1  Presentation  of  Data  and  Results 

Data  presented  here  will  demonstrate  that  the  wave  interaction  experiment 
has  sufficient  signal  to  noise  ratio  to  detect  dynamic  perturbations  in  the 
D-region.  Figure  2. 14  shows  one  hour  of  data  taken  August  2,  1976,  at  around 
noon  time,  at  82  km,  and  one  hour  of  noise  from  the  same  recorder  channel  taken 
at  a time  when  ionospheric  conditions  were  almost  the  same  as  when  the  data 
were  taken.  Conditions  were  very  good  in  both  cases,  based  upon  operator 
judgement  regarding  the  stability  of  the  echo  and  the  absence  of  interfering 
signals . 
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The  fluctuations  in  the  noise  curve  are  not  due  to  system  noise,  but  rather 
to  ionospheric  noise  and  echo  fluctuations  as  described  earlier.  System  noise 
is  about  another  factor  of  five  below  the  noise  shown  here  and  is  thus  totally 
insignificant.  It  is  clear  from  Figure  2.14  that  the  fluctuations  on  the  data 
are  much  greater  than  what  may  be  explained  by  noise.  Figure  2.15  shows  the 
same  comparison  for  data  from  56  km,  where  the  data  is  much  smaller.  The  noise 
from  that  channel  is  approximately  the  same  as  in  all  channels,  so  the  signal 
to  noise  ratio  is  not  as  good  at  low  heights.  Nonetheless,  it  is  still  adequate. 

The  power  spectra  for  the  two  heights,  data  and  noise,  are  shown  in  Figures 
2.16  and  2.17.  Since  the  power  spectrum  is  proportional  to  amplitude  squared,  the 
signal  to  noise  ratio  is  exaggerated  and  when  the  two  are  plotted  on  the  same 
scale,  the  noise  spectrum  cannot  be  distinguished  from  the  base  line  at  82  km, 
and  is  still  quite  small  at  56  km.  Conditions  could  deteriorate  considerably 
before  noise  becomes  a problem.  Under  certain  unusual  conditions,  noise  can 
be  as  much  as  40  db  greater  in  effect  than  shown  here.  Such  conditions  occur 
when  the  echo  is  both  weak  and  unstable,  and  background  noise  is  high.  Such 
conditions  are  unusual,  however,  and  the  experiment  usually  operates  with  a 
satisfactory  signal  to  noise  ratio. 

Power  spectra  cannot  prove  the  existence  of  waves,  but  those  presented  here 
suggest  there  may  be  periods  present  which  are  significantly  above  the  background 
level.  It  is  possible  that  such  periods  are  the  results  of  waves,  although 
gravity  waves  are  expected  on  the  average,  to  have  a continuous  spectrum  without 
significant  peaks.  However,  any  specific  sample  could  exhibit  such  peaks.  It 
is  interesting  that  the  noise  spectrum  at  56  km  shows  a few  significant  peaks 
also.  This  can  be  explained  as  some  sort  of  periodic  behavior  in  the  amplitude 
of  the  wanted  echo.  It  is  not  known  whether  or  not  the  cause  is  wavelike 
behavior  in  the  E-region  at  the  reflection  height. 


Figure  2.15  Comparison  of  data  and  noise  at  56  km 
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To  illustrate  the  capabilities  of  the  entire  system  several  sets  of  data 
and  spectra  are  shown  in  Figures  2.18  - 2.23.  These  represent  data  from  what 
are  known  as  all-day  runs,  since  the  period  of  data  collection  covers  most  or 
all  of  the  daylight  hours.  The  data  and  spectra  curves  indicate  the  possibility 
of  wavelike  behavior.  The  next  subsection  discusses  the  steps  necessary  for  a 
more  detailed  analysis  of  this  data. 


2.3.2  A Summary  of  Present  Work 

As  discussed  in  2.1  of  this  report,  data  taken  at  a particular  tu  is  not  a 
result  of  ionospheric  parameters  from  that  height  alone.  The  range  of  influence 
extends  from  that  height  down  to  the  base  of  the  ionosphere,  hut  a detailed  study 
of  the  experiment  shows  that  heights  farther  from  h^  exert  less  influence.  This 
is  a result  of  the  very  rapid  cooling  times  in  the  lower  D-region.  The  range  of 
influence  also  extends  upward  from  h^  due  to  the  finite  widths  of  the  pulses. 

None  of  this  causes  height  resolution  to  be  completely  lost,  but  of  course,  it 
is  not  possible  to  fully  recover  the  one  kilometer  resolution  which  is  the 
absolute  limit  due  to  data  spacing. 

The  technique  for  height  resolution  recovery,  which  has  been  developed  and 
now  being  programmed  is  basically  a deconvolution  of  finite  wanted  and  disturbing 
pulse  data  into  that  of  very  narrow  wanted  and  disturbing  pulses.  For  such  an 
approach,  it  is  necessary  to  extend  the  data  profile  below  the  lowest  altitude 
at  which  data  has  been  taken  by  assuming  that  the  data  tapers  rapidly  to  zero. 
This  extrapolation  agrees  quite  well  with  the  actual  physical  situation  because 
the  increasing  collision  frequency  at  lower  altitudes  does  cause  the  data  values 
to  decrease  rapidly  with  height.  It  may  be  shown  that  the  wave  interaction 
coefficient  at  a designated  height  of  a finite  wanted  pulse  is  the  mean  of  the 
sum  of  the  coefficients  of  a set  of  very  narrow  pulses  with  height  range  covering 
the  finite  pulse.  The  deconvolution  process  has  to  be  applied  to  both  the  wanted 
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Power  spectra  of  three-dimensional  data  (8121 


Figure  2.22  Three-dimensional  wave  interaction  data  (electron  precipitation  event 
March  27,  1976,  State  College) 
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Figure  2.23  Power  spectra  of  three-dimensional  data  (electron  precipitation 
event,  March  27,  1976,  State  College) 
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and  the  disturbing  pulses.  In  the  case  of  the  wanted  pulse  if  the  finite 
pulse  coefficients  are  known  every  kilometer,  the  coefficients  for  the  narrow 
pulses  may  be  found  by  starting  at  the  base  of  the  ionosphere  where  the 
coefficients  are  assumed  to  be  zero  and  working  upward. 

The  finite  D pulse  may  be  treated  in  a similar  manner,  and  using  a similar 
model  tapering  the  data  points  to  zero,  the  coefficients  for  narrow  wanted  and 
disturbing  pulses  may  be  recovered  from  those  for  narrow  wanted  and  wide 
disturbing  pulses.  The  accuracy  of  these  two  processes  is  limited  by  the 
accuracy  of  the  assumed  model  at  low  heights,  and  also  by  the  accuracy  of  the 
data  points. 

Now  that  the  data  has  been  processed  into  the  form  where  it  may  be  modeled 
by  the  conceptually  simple  narrow  wanted,  narrow  disturbing  pulse  experiment,  it 
is  fairly  simple  to  find  a process  for  deriving  quantities  which  are  proportional 
to  parameters  at  one  height  only.  Letting  the  wave  interaction  coefficient  for 
narrow  pulses,  at  h^  be  given  by  T(h^),  it  may  be  shown  that 

hi 

T(h  ) = / W(h)  J(hi,h)  dh  2-27 

o 

where  W(h)  is  a quantity  which  depends  only  upon  ionospheric  parameters  at  h 

assuming  that  the  power  absorbed  in  the  lower  D-region  from  the  D pulse  transmission 

-2  Gv(hi-h) 

is  negligible,  and  J(h.,  h)  = exp  ( — ) where  v is  the  electron-neutral 

1 c 

collision  frequency  and  G is  the  energy  loss  factor. 

The  above  integral  may  be  approximated  by  a summation  with  one  kilometer  steps. 
Making  the  same  assumption  as  before  with  the  data  tapering  to  zero  the  quantity 
W(h)  may  be  found  for  a given  assumed  height  profile  of  Gv.  Note  that  each 
T(h^)  is  a sum  of  W(h) 's  weighted  by  the  appropriate  J's.  Beginning  at  low 
altitudes  where  the  data  values  are  zero  and  working  up  will  yield  a set  of  W's 
whose  accuracy  depends  upon  the  accuracy  of  the  assumed  Gv  profile.  Since  the 


decay  is  rather  rapid  at  low  altitudes,  some  inaccuracy  in  the  Gv  profile 
can  be  tolerated  while  still  achieving  reasonably  good  height  resolution. 


The  next  step  is  to  normalize  the  data  so  that  changes  can  be  expressed 
as  fractional  variations.  Finally,  cross-correlation  and  co-spectrum  analyses 
may  be  performed. 


, t ssk*  ss+m.  * . . ; nth* 


% 


53 


3.  Intercomparison  of  Electron  Density  Synthesis  Methods 

The  difficulties  of  obtaining  accurate  electron  density  height  profiles 
from  groundbased  experimental  data  profiles  arise  from  the  facts  that  the 
data  profiles  are  usually  an  integral  effect  over  a wide  height  range  and  the 
data  are  contaminated  with  noise.  Therefore  it  is  important  to  somehow 
enhance  the  reliability  of  synthesized  electron  density  profiles  through 
statistical  means.  The  two  methods,  Dale's  and  Shellman's,  to  be  compared 
in  this  section  both  have  common  aim  of  minimizing  the  afore-mentioned 
difficulties.  As  will  be  seen  from  the  detailed  discussions  given  in  the 
following  subsections.  Dale's  method  achieves  this  aim  by  applying  overdeter- 
ministic  approach  in  which  the  data  at  larger  number  of  heights  are  used 
compared  to  the  number  of  heights  where  unknown  electron  densities  are  being 
determined.  On  the  other  hand  Shellman's  method  achieves  the  same  aim  by 
balancing  the  trade-off  between  the  smoothness  of  electron  density  profiles 
and  the  fit  to  the  data  as  described  later  in  this  section. 

Throughout  this  section  the  amplitude  interaction  coefficient  T^  is 
used  as  an  example  and  both  methods  assume  all  parameters,  such  as  collision 
frequency,  electron  temperature,  energy  loss  coefficient,  etc.,  as  known; 
thus  treating  electron  densities  as  the  only  unknown.  In  subsection  3.1  and 
3.2,  theoretical  bases  for  both  methods  are  presented;  whereas  in  3.3,  test 
results  from  the  two  methods  are  intercompared. 

3.1  Dale  Synthesis  Method 

Since  the  basic  theory  of  the  wave  interaction  experiment  involved  in 
this  intercomparison  has  been  well  documented  in  previous  papers,  only  the 
expressions  appropriate  for  introducing  synthesis  method  will  be  cited  here. 
Detailed  development  of  the  theory  can  be  found  in  Lee  and  Ferraro  (1969) 


and  Dale  (1973). 


As  mentioned  earlier,  the  so-called  over-deterministic  synthesis 


method  is  basically  a method  of  inverting  a data  height  profile  to  yield 

an  electron  density  height  profile  using  larger  number  of  data  points  than 

the  number  of  heights  where  electron  densities  are  determined. 

The  amplitude  and  phase  interaction  effects,  denoted  from  here  on  by 

T.  and  T , respectivelv,  are  simply  expressed  in  terms  of  the  perturbation  of 
A <p 

the  imaginary  (x)  and  real  (g)  parts  of  the  Sen-Wyller  refractive  index  for 
the  wanted  signal  as 


T 


h. 


= j J1  A n (h)  S dh 


= 1 


/ s 


3. W 
3v 


3v 


Av  dh 


or 


X = / 6 v — 

A 1 D 3v  v 
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r a 3u  av  ,, 

* 6 37  V V dh 
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Av 


dh 


3-1 


3-2 


where  R'  is  the  reflection  coefficient  of  the  disturbed  sequence  of  wanted 

pulses  and  R is  that  of  the  undisturbed  sequence  of  wanted  pulses.  £ is  the 

free  space  propagation  constant  for  wanted  pulses  and  v is  the  most  probable 

collision  frequency  of  Sen-Wyller  formula  for  index  of  refraction.  A 

different  definition  of  T,  for  example  T = (R'  - R)/R,  will  result  in  an 

opposite  sign  in  the  expressions  of  T and  T . However,  as  long  as  the  polarity 
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of  the  quantities  T and  T measured  by  the  detecting  system  is  made 
A 

consistent  with  the  definition,  use  of  either  definition  is  acceptable. 

It  should  be  noted  that  the  quantity  can  be  either  negative  or  positive 

dV 

depending  on  height  range.  It  is  typically  in  the  case  of  the  ordinary 
mode  wanted  wave,  negative  below  around  65  km  and  positive  above,  thus 
creating  a crossover  height  between  the  negative  and  positive  amplitude 
interaction  data.  Both  y and  x involve  the  electron  density,  Ng,  and  the 
collision  frequency,  v;  and,  in  summary,  Equations  3-1  and  3-2  can  be 
expressed  in  the  form  of 


XA,  <V  = 


(h,hi>  F(h)  dh 


3-3 


where  F(h)  is  the  only  term  which  contains  the  unknown  electron  density 
N . h is  the  base  of  the  ionosphere  and  h.  is  the  height  of  interaction 

6 D 1 

where  the  upgoing  disturbing  wave  meets  the  downgoing  wanted  wave  in  the 

wave  interaction  experiment.  Other  variables  of  height  such  as  collision 

frequency,  which  is  contained  in  0.  , (h,h.)  and  F(h)  are  considered  to  be 

A,<t>  i 

known  in  this  synthesis  method.  Subscripts  A and  <fc  denote  amplitude  and 
phase  interaction  respectively. 

F(h)  is  the  portion  of  the  integrand  that  contains  all  the  electron 
density  terms  in  the  form  of 


F(h) 


exp 


u 


f(z)  dz 


B 


3-4 


where  z is  a dummy  variable  for  height. 

The  general  approach  of  the  synthesis  is  to  solve  for  F(h)  using  the 
system  of  simultaneous  equations  obtained  by  forming  Equation  3-3  for  each 
data  point  TCIk)  where  the  number  of  TOk)  is  larger  than  that  of  unknown  F(h). 
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The  first  step  of  the  technique  is  to  solve  for  F(h).  Since  linear 
interpolation  is  used  on  F,  as  a function  of  height  h,  below  to  over- 
determine the  system,  it  is  desirable  that  the  unknown  function  F(h) 
be  fairly  linear  over  the  interpolated  range.  To  accomplish  this,  an 
estimate  of  the  electron  density  profile  is  made  and  substituted  into 
3-4.  The  resulting  known  function,  referred  to  as  F9(h),  is  multiplied 

by  an  unknown  correction  function,  F .(h) , and  the  product  is  made  equal 

*■  « 

to  F(h).  That  is, 

F (h)  = F^h)  F2(h)  3-5 


This  is  substituted  into  3-3.  If  the  estimate  is  close  to  the  actual 
profile  the  F(h)  will  be  approximately  linear  over  small  height  increments. 
Numerous  tests  indicated  that  such  an  indirect  approach  as  using  an  initial 
estimate  of  an  electron  density  profile  improved  the  results  over  those 
found  by  solving  for  F(h)  directly.  Each  of  these  F's  are  assumed  constant 
over  horizontal  slabs  of  ionosphere.  This  reforms  Equation  3-3  as 


T(I)  = l F (K)  F (K)  / 

K 


Kth  slab 


G(h,  hi>  dh 


3-6 


It  was  found  that  integrating  G(h,h^)  over  the  slab  gave  much  better 
results  than  assuming  it  to  be  constant  over  the  slab. 

Now,  define  the  elements  of  a matrix  (A°)  to  be 


A°  (I,  K)  = F,(K)  / 


Kth  slab 


GCh.hJ  dh 


3-7 


where  (A  ) has  nonzero  elements  in  positions  consistent  with  the  range  of 
the  sum  in  Equation  3-6  or  the  limits  of  the  integral  in  Equation  3-3. 
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The  system  can  now  be  expressed  in  matrix  form  as 


UJ 


La°: 


3-8 


where 

r 


lTJ 


(n  x 1)  matrix  with  n representing  the  number  of  data  points 


r - 


jF 

L J 


(m  x 1)  matrix  with  m representing  the  number  of  heights  where 


N is  determined 
e 


rA°' 


(n  x m)  matrix 


In  order  to  introduce  over-deterministic  feature,  let  be  formed 


by  the  linear  interpolation  of  the  elements  of  R , where  the  ratio 


of  the  number  of  elements  of  to  those  of  R is  the  degree  of 


over-determination.  This  can  be  expressed  as 


>J 


Is'  r-~ 


IS  1 R 

— j 


3-9 


The  system  of  equations  to  be  solved  is  then 


!T 


[a°  rs 


3-10 


or,  multiplying  Fa°J  and  fsj 


Tj  = jAS1  rR 


3-11 


Lastly,  to  eliminate  possible  negative  values  for  the  elements  of  jR  , 
define  a new  vector  R'  such  that 


R(J)  « exp  { R ’ (J)  } 
where  J is  an  index. 


3-12 
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The  system  of  equations  defined  by  Equations  3-11  and  3-12  are  now 

solved  by  a nonlinear,  least  square  algorithm  of  Marquardt  (1963,  1971),  the 

r ' r - 

result  of  which  is  the  vector  R'  . R is  formed  using  Equation  3-12 

and  in  turn  [ F^  by  Equation  3-9  and  F by  Equation  3-5. 

All  that  remains  is  to  solve  for  , knowing  F and  using  Equation  3-4 
which  is  rewritten  in  the  form 

j / h ■)  *1 

N = F exp  ^ - / N f (z)  dz  - 3-13 

e i i i e 

L \ hB  J 

This  equation  is  solved  using  an  induction  method  starting  at  the  lowest 
slab  and  solving  for  the  electron  density  at  each  skab  using  Wegstein 
iteration  scheme  (Lance,  1960  and  Wegstein,  1960)  which  solves  equations 
of  the  form  x = f (x) . 

The  effect  of  a finite  width  wanted  pulse  in  the  wave  interaction 
experiment  can  easily  be  included  by  using  the  convolution  integral: 


T (x)  = — / T(x')  dx' 

w CT  /n 

ct/2 

where  x is  the  width  of  the  wanted  pulse  or  receiver  gate  in  seconds.  This 
causes  only  a modification  to  the  matrix  Ha°  and  does  not  alter  the  basic 
technique.  In  this  intercomparison,  the  effects  of  a finite  width  wanted 
pulse  are  fully  taken  into  account. 

The  theory  of  synthesis  method  presented  here  is  equally  applicable  to 

amplitude  and  phase  interaction  data  T^  and  T^ . However,  the  intercomparison 

tests  were  conducted  concentrated  on  T.  onlv. 
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3.2  Shellman  Synthesis  Method 

Since  basic  wave  interaction  theory  has  been  described  in  summary 
form  in  3.1  along  with  the  Dale  synthesis  method,  only  the  key  factors 
depicting  the  basic  concept  of  the  Shellman  Synthesis  Method  will  be 
summarized  here.  The  details  of  the  method  could  be  found  in  Shellman 
(1973). 

In  this  method,  in  order  to  parameterize  the  electron  density  profile, 
it  is  represented  by  a series  of  short  exponential  segments.  The  common 
logs  of  the  electron  densities  at  the  heights,  hj  , where  the  segments 
join,  are  taken  to  be  the  unknown  parameters  of  the  profile  and  are 
denoted  by  = log  Ne(h  ),  where  j = 1,  2,  . . .,n. 

Basic  approach  in  minimizing  experimental  errors  is  to  search  for  the 
smoothest  possible  electron  density  profile  which  provides  exceptable  fit 
to  the  experimentally  measured  data.  The  quantitative  scheme  of  simultan- 
eously satisfying  the  smoothness  and  the  data  fit  is  to  minimize  the 
quantity 

f = c + Xs  3-14 


where 


a2  2 

- / ( V ) 

dz 


3-15 


is  the  measure  of  curvature  in  the  electron  density  profile  with 

a = log  N and 

e 


m / t(h.)  - T(h . ) j 2 
= l i — . 
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i=l 
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is  the  measure  of  deviation  of  data  value,  t(h^),  computed  from  the 
assumed  electron  density  profile,  from  experimentally  measured  data 
T(h^)  with  as  standard  deviation.  The  equation  3-14  could  take  the 
form  of 


f = c + 4 X s 


to  simplify  mathematical  manipulations  or  the  form  of 


2 

f = c + d + % X s 


in  order  to  include  a damping  term  d where 


d = K , / "a (log  N )]  2 dz 
d ' l e -i 


where  is  some  chosen  constant  and  A (log  N ) is  the  change  in  the  profile 
at  height  z corresponding  to  the  increment  AX. 

General  steps  taken  in  the  search  scheme  is  first  to  find  the  electron 
density  profile  corresponding  to  X = 0,  which  is  the  best-fit  exponential 
profile;  and  then  find  profiles  for  successively  larger  values  of  X as 
details  are  added  to  the  profile  with  the  progress  of  search.  The 


increments  in  X depend  on  the  linearity  of  the  t(10  with  respect  to  . 

For  the  present  intercomparison  purpose,  equation  3-14  was  used  while 
the  damping  term  d is  included  only  in  the  process  of  determining  the  best- 


fit  exponential  profile.  The  advantage  in  starting  with  the  best-fit 

StOiJ 

exponential  profile  is  that  the  derivatives  — — are  better  behaved 

for  profiles  of  simpler  form.  This  is  an  important  factor,  as  can  be 
seen  from  the  search  scheme  described  below. 
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The  first  step  in  obtaining  a solution  to  Equation  3-14  is  to 
differentiate  this  equation  with  respect  to  a_.'s,  yielding 


^c_  + x Js_ 

3a . 3a . 

3 3 


» 9t(h  ) 

taT  + 2.i,  “to“  "°1 
3 I 1=1  3 


■ [t(h.)  - T(h.)]  = 0 

- 


This  expression,  however,  cannot  be  solved  directly  since  t(h.)'s  are  not 

3t(h.)1 

known  unless  a.'s  are  known  and  since  the  derivatives  — ; are  functions 

j da . 

3 

of  a^  . Nevertheless,  solution  equations  can  be  given  for  the  a.'s 

in  terms  of  the  parameters  a.  = a of  some  assumed  starting  profile,  and 

3 J 3t(h  ) 

the  computed  values  t(h. ) = tn(h  ) and  derivatives  — , for  the  assumed 

1 i dot . 

3 

profile.  These  solution  equations  can  be  used  if  the  starting  profile 
and  solution  profiles  are  sufficiently  similar. 

In  application  the  curvature  function  c of  Equation  3-15  is  expressed 
in  terms  of  summation  over  profile  segments.  This  becomes 


a . . - a . a . - a . , . 

3-1  3 3 3+1 


n-1  h - h . h . 

y J-1 j l 

A _ <4h>j 


m -i 

~.hi.+ i 


where  the  expression  in  brackets  represents  a second  derivative  at  height 

h . and 
3 


(At,).  - (hj.j  - h.+1>/2 


is  the  height  increment  over  which  the  j-th  second  derivative  is  considered 
to  apply.  The  slab  summation  form  of  c approaches  the  integral  of  Equation 
3-15  as  the  slab  is  made  sufficientlv  thin. 
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The  search  scheme  described  above  has  a tendancy  to  give  successively 
smaller  value  of  s,  in  the  process  of  iteration,  at  a substantial  rate  until 
it  reaches  a threshold  beyond  which  the  rate  of  reduction  in  s become 
insignificant.  The  filial  electron  density  profile  is  chosen  at  this 
threshold  point. 

3.3  Results  of  Comparison 

Intercomparison  of  the  two  electron  density  synthesis  method  was 
initially  made  with  the  data  generated  from  known  electron  density  profiles. 
Then  as  a final  test  the  two  methods  were  applied  to  real  data  obtained 
from  the  wave  interaction  experiment  conducted  at  the  Ionosphere  Research 
Laboratory,  PSU. 

Eighteen  samples  obtained  over  a five-hour  period  on  January  31,  1976 
were  synthesized  by  both  methods.  Figures  3.1  through  3.3  show  three 
typical  comparisons,  while  fifteen  others  are  presented  as  Figures  1-A 
through  15-A  in  the  appendix  for  reference. 

As  can  be  seen  from  the  figures,  the  results  of  two  different  methods 
are  in  excellent  agreement.  However,  the  Dale  synthesis  method  seem  to  give 
somewhat  more  detail  in  electron  density  profiles.  This  could  be  due  to 
the  fact  that  the  Shellman  Synthesis  Method  has  a built-in  feature  of  balancing 
the  trade-off  between  the  smoothness  of  profiles  and  the  data  fit  which  tends 
to  give  smoother  profiles. 
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Figure  3.2  Intercomparison  of  synthesized  electron  densities 
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4.  Routine  Operation  for  Monitoring  D-Region  Temporal  Variability 

As  stated  earlier  in  the  introduction,  this  section  describes  the 
capability  of  IRL's  digitalized  sophisticated  wave  interaction  experimental 
facilities  in  monitoring  and  studying  temporal  variations  in  the  D-region 
ionosphere . 

The  wave  interaction  data,  which  are  digitally  recorded  on  magnetic 
tapes  covering  28  D-region  heights  simultaneously  once  every  300  milli- 
seconds, can  be  processed  in  a diversified  manner.  It  is  best  to  describe 
the  processes  center  around  the  flow  chart  shown  in  Figure  4.1.  At  present 
the  Penn  State  Facility  is  not  directly  linked  to  a computer  although  in 

i 

principle  it  can;  it  has  not  been  necessary  for  us  to  do  so  to  meet  our 
scientific  objectives.  We  do  still  have  rather  fast  "tum-around"  with  the 
procedure  given  in  Figure  4.1.  A wave  interaction  data  tape  is  quickly 
processed  at  Stage  I to  give  a "quick  look"  at  data  profiles;  this  takes 
only  several  minutes  after  delivering  the  tapes  to  Computation  Center. 

Some  averaging  is  done  at  this  stage  since  during  an  8 hour  run,  for 
example,  there  are  almost  2.7  x 10^  data  values  for  the  height  range  50  to 
90  km.  Stage  I then  gives  us  an  opportunity  to  examine  the  quality  of  data 
and  from  our  own  experience  allows  us  to  specify  the  state  of  the  D-region 
ionization.  A condensing  program  allows  us  average  data  profiles  corres- 
ponding to  time  periods  of  several  seconds  to  several  minutes  depending 
upon  what  D-region  feature  is  being  investigated.  Minutes  would  be  sufficient 
for  evaluating  diurnal  features  or  relativistic  electron  precipitation  events 
whereas  seconds  would  be  needed  for  X-ray  flare  events,  for  instance.  Stage  II 
conveniently  displays  the  data  in  "3-dimensional"  format  as  in  Figure  4.2 
for  the  quiet  D-region  of  January  31,  1976.  It  is  evident  from  Figure  4.2 
that  there  is  a diurnal  effect  as  well  as  irregular  behavior  in  terms  of 
"wave  structure"  present.  Stage  II  output  takes  1-3  hours  to  obtain  dependent 
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upon  computer  availability.  If  atmospheric  wave  structure  information  is 
required,  then  the  output  of  the  condensing  program  is  put  through  a fast 
Fourier  transform  into  Stage  III  for  a display  of  the  spectral  data. 

Figure  4.3  is  such  a "3-D"  presentation  corresponding  to  the  data  of  Figure 
4.2.  Note,  periods  of  atmospheric  waves  range  from  7 hours  to  .35  hours  with 
correlation  of  the  peaks  with  altitude.  An  output  from  the  condensing 
program  can  be  transferred  to  the  Dale-Shellman  synthesis  routine  described 
earlier  in  Section  3 for  generation  of  electron  density  profiles  versus 
time.  It  should  be  noted  that  all  the  results  presented  in  this  section 
were  obtained  from  Dale's  synthesis  method.  All  of  these  stages,  of  course, 
contribute  to  the  immense  flexibility  of  the  experiment  to  support  a number 
of  programs  suggested  in  Stage  VI. 

Figure  4.4  shows  the  diurnal  variation  of  electron  densities  at  various 
altitudes  for  the  quiet  day  of  January  31,  1976. 

Figure  4.5  shows  a 3-D  display  of  all-day  wave  interaction  data  on 
March  27,  1976  during  the  recovery  phase  of  a geomagnetic  storm  which  followed 
auroral  activity.  The  important  features  are  the  extremely  increased 
negative  cross-modulation  at  low  heights  during  the  earlier  times  which 
indicate  production  is  not  solar  and  probably  electron  precipitation;  lack 
of  a clear  diurnal  feature  as  in  Figure  4.2  since  electron  precipitation 
dominated  over  solar  control.  Figure  4.6  displays  a few  electron  density 
profiles  which  correspond  to  larger  densities  of  such  disturbed  day 
compared  to  these  of  normal  quiet  conditions. 

It  is  evident  from  the  illustrations  given  in  this  section  that  the 
all-digitalized  wave  interaction  experimental  facilities  now  possess  the 
time  and  height  resolutions  capable  of  determining  temporal  variations  of 
electron  densities  as  well  as  studying  wave  motions  in  the  D-region  iono- 
sphere under  normal  quiet  and  disturbed  conditions. 
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Figure  4.3  Power  spectra  of  three-dimensional  data  (quiet  day 
January  31,  1976,  State  College) 
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This  appendix  contains  fifteen  additional  figures  showing  the 
results  of  intercomparison  between  the  Dale  and  the  Shellman  electron 
density  synthesis  methods  described  earlier  in  Section  3. 
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